clear;
clc;

cd('specify directory here');
load('EdEx_matlab_years');

global Y T1 T2 SM SR Draws M ll 
options=optimoptions('fminunc','Algorithm','quasi-newton','Display','iter','PlotFcn',{@optimplotx,@optimplotfval,@optimplotfirstorderopt});

rng(10010);
Drawso=normrnd(0,1,6063,500);
M=500;

SM=(SM-mean(SM))/std(SM);
SR=(SR-mean(SR))/std(SR);

Y=Yo; T1=T1o; T2=T2o; SM=SMo; SR=SRo; Draws=Drawso; 

theta2=[0;0;1;1;1;1;0;1;1;0.6;0.6;0;0;1;1;0.2;0.2];
eps=1e-5;
no = 0.1;
while no>eps
    theta = theta2;
    [theta2, FVAL]=fminunc(@(param)LogLikelihood_years(param),theta,options);
    no=norm(theta2-theta);
end

PARAM = theta2;
var = ERRORS_years(PARAM);
se = sqrt(diag(var));
